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In a quantizing magnetic field, the two-dimensional electron (2DEG) gas has a rich phase diagram 
with broken translational symmetry phases such as Wigner, bubble, and stripe crystals. In this 
paper, we derive a method to get the dynamical matrix of these crystals from a calculation of the 
density response function performed in the Generalized Random Phase Approximation (GRPA). 
We discuss the validity of our method by comparing the dynamical matrix calculated from the 
GRPA with that obtained from standard elasticity theory with the elastic coefficients obtained from 
a calculation of the deformation energy of the crystal. 



PACS numbers: 73.20.Qt,73.21.-b,73.43.-f 



I. INTRODUCTION 

Theoretical calculations show that, in the presence of 
a perpendicular magnetic field, a two-dimensional elec- 
tron gas (2DEG) should crystallize below a filling factor 
v ^ 1/6. 5i. Several experimental groups have reported 
transport measurements indicative of this electron crys- 
tallization when the filling factor of the lowest Landau 
level is decreased below v = 1/5. These measurements 
include the observation of a strong increase in the diago- 
nal resistivity p^i-^, non-linear I — V characteristics, and 
broadband noise. All these observations have been in- 
terpreted as the pinning and sliding of a Wigner crystal 
(WC)-. Moreover, microwave absorption experiments'^ 
have also detected a resonance in the real part of the lon- 
gitudinal conductivity, cfxx (i-^), that has been attributed 
to the pinning mode of a disordered Wigner crystal. The 
vanishing of the pinning mode resonance at some critical 
temperature Tm (v) has been used to derive the phase 
diagram of the crystal^ in the quantum regime where 
the kinetic energy is frozen by the quantizing magnetic 
field. Similar microwave absorption experiments also 
showed a pinning resonance at higher filling factors close 
to v = 1,2,3 where the formation of a Wigner solid is 
expected in very clean samples"'^*?. Finally, in Landau 
levels of index > 1, a study of the evolution of the pin- 
ning mode with filling factor reveals several transitions of 
the 2DEG ground state from a Wigner crystal at low v to 
bubble crystals with increasing number of electrons per 
lattice site as v is increased, and into a modulated stripe 
state (or anisotropic Wigner crystal) near half filling^. 

In earlier works^iii^, some of us have studied several 
crystalline states of the 2DEG using a combination of 
Hartree-Fock (HFA) and generalized random-phase ap- 
proximations (GRPA). In these works, the energy and 
order parameters of the crystal were calculated in a self- 
consistent HFA while the collective excitations were de- 
rived from the poles of density response functions com- 
puted in the GRPA. This microscopic approach (HFA 
-I- GRPA) works well at zero temperature but is diffi- 
cult to generalize to consider finite temperature effects 



or to include quantum fluctuations beyond the GRPA. 
Finite-temperature or quantum fluctuations effects (not 
already included in the GRPA), are most easily com- 
puted by writing down an elastic action for the system. 
For a crystaline solid, this requires the knowledge of the 
dynamical matrix (DM) or, equivalently, of the elastic 
coefficients of the solid. 

A direct way to obtain these elastic coefficients is to 
compute the energy required for various static defor- 
mations of the crystal. Using elasticity theory, each 
deformation energy AEi can be written in the form 
AEi — jCiU^ where is a parameter characterizing 
the amplitude of the deformation and Ci is generally a 
combination of elastic constants. In the limit 0, one 

can obtain the elastic coefficients by computing the de- 
formation energy of one or more static deformations and 
using the known symmetry relations between the elas- 
tic constants. Alternatively, one can obtain a DM from 
the GRPA density response function much more directly 
without the need to compute the elastic coefficientai^. 

In this paper, we compare the DM obtained from these 
two methods (deformation energy and GRPA) in order 
to find the range of validity as well as the limitations of 
the GPRA approach. We first consider the simple case of 
an isotropic (triangular) Wigner crystal before tackling 
the more complex anisotropic Wigner crystal^- or stripe 
phase that occurs near half-filling in the higher Landau 
levels. We show that although the GRPA method gives 
a good description of the qualitative behavior of the DM 
as a function of filling factor, its quantitative predictions 
must be used with caution. As we show below, an averag- 
ing procedure must be applied to the method in order to 
obtain a DM in the GRPA that compares favorably with 
the one obtained by computing the deformation energy. 

Our paper is organized as follows. In Section II, we de- 
fine the elastic constants needed to build an elastic model 
for the Wigner and stripe crystals. We then explain, in 
Section HI, how these elastic constants can be derived by 
computing the deformation energy of the crystals in the 
HFA. In section IV, we summarize the GRPA method 
of obtaining the dynamical matrix. Our numerical re- 
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suits for the WC are discussed in Section V and those for 
the stripe crystal in Section VI. Section VI contains our 
conclusions. 



II. ELASTIC CONSTANTS AND DYNAMICAL 
MATRIX 



We describe the elastic deformation of a crystal state 
by a displacement field u (R) defined on each lattice site 
R. The Fourier transform of this operator is given by: 



u(k) 




ik-R 



u(R), 



(1) 



where is the number of lattice sites. In two dimen- 
sions, the general expression for the deformation energy 
of a crystal requires the use of 6 elastic coefficients Cij 
and is given, in the continuum limit, by the following 



expressioi 



AE = - 
2 



i J dr[ciiel,^ + 4:CGQely + AcQ2e^^yeyy] ( 
- / dr [2ci2 + 4ci6 



where 



If dua (r) 



2 V ar, 
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dup (r) 



(3) 



is the symmetric strain tensor. 

The Wigner and bubble crystals have a triangular lat- 
tice structure for which the following equation holds: 



Cll = C22 = 2C66 + Ci2. 



(4) 



For such a triangular structure, the elastic energy in the 
long- wavelength limit can be written in a form that con- 
tains only two elastic coefficients, namely: 

= \ J dr [ci2 {el, + ely + 2e,^,ey,y) 

+ 2c66 {el, + ely + 2ely)] . (5) 

The anisotropic stripe state can be seen either as a 
centered rectangular lattice with two electrons per unit 
cell or as a rhombic lattice with one electron per unit 
cell with reflection symmetry in both x and y axis. The 
deformation energy is given by 

AS = ^ y rfr [cnel, + 'Iceeely ] 

+ - dr [2ci2e + C22eyy] ■ 

In this paper, we assume that the stripes are aligned 
along the y axis. 

The above formulation of elasticity theory assumes 
short-range forces only. For the electronic crystals that 



we consider, these forces are of coulombic origin i.e. the 
hamiltonian of the crystal contains only the Coulomb in- 
teraction between electrons and the kinetic energy which 
is frozen by the quantizing magnetic field. Both the di- 
rect (Hartree) and exchange (Fock) terms are considered 
by the Hartree-Fock approximation as we explain in the 
next section. To take into account in the elasticity theory 
the long-range part of the Coulomb interaction present 
in a crystal of electrons, it is necessary to add to /S.E the 
deformation energy lS.Ec given by 



AEr 



dv 



,6n (r) Sn (r') 



(7) 



7re^ ^ Sn (q) Sn (— q) 



S ^ 
q 



where S is the area of the crystal, Sn (q) = 
/ dxe^'^'^'^Sn (r) is the Fourier transform of the change 
in the electronic density and n, is the dielectric constant 
of the host semiconductor. We consider the positive 
background of ionized donors as homogeneous and in- 
ert so that no linear term in Sn (r) is introduced by the 
Coulomb interaction. 

To define a dynamical matrix, we assume that the crys- 
tal can be viewed as a lattice of electrons with static form 
factor h (r) on each crystal site (with the normalisation 
/ drh{r) = 1). The time-dependent density can then be 
written as 



,(r,t) =^/i(r-R-u(R,i)), 



(8) 



R 



and, to first order in the displacement field, we have for 
a density fluctuation 



Sn(k + G,t) = -ih{k 



G).u(k,t), (9) 



where G is a reciprocal lattice vector and k a vector in 
the first Brillouin zone of the crystal. It follows that we 
can write the Coulomb energy as 



AEc = TTTioe 



q 



|/t(q)q-u(q) 

Kq 



(10) 



where no = Ns/ S is the average electronic density. 
We pause at this point to remark that the form factor 



hiq) 



(with £ — ^ hc/eB the magnetic length. 



B being the applied magnetic field) in Eq. (fTO|) renders 
the summation over the wavectors rapidly convergent. 
Our Hartree-Fock calculation of the ground-state energy 
of the electronic crystals as well as our CRPA calculation 
of the dynamical matrix also involve summations over 
reciprocal lattice vectors G of some functions weighted 
by h (G). If the magnetic field is not too strong, we can 
perform these summations directly. There is no need to 
use Ewald's summation technique as is the case if one 
works with a crystal of point electrons. Of course, as 
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the filling factor p —f 0, the magnetic length i ^ so 
the electrons behave more and more like point particles 
and the convergence is lost. In all cases that we consider, 
the summations involved are rapidly convergent because 
we restrict ourselves to filling factors i' — 27rno^^ ^0-1 
where £/ao is sufficiently large for ^ to be small 
(ao being the lattice constant). The cutoff in G is choosen 
so that the summations are evaluated with the required 
degree of accuracy. 

The total deformation energy, which we now write as 
A.Ex, now includes the long-range Coulomb interaction 
and can be written in the form: 

A^^T = ^ ^ "a (k) (k) (-k) , (11) 

k 

where we have introduced the dynamical matrix: 

For the triangular lattice, a comparison of Eqs. pT|) and 
([5|) gives the dynamical matrix (to order k'^) as: 

D^^^ (k) = 1 [{ci2 (k) + 2c66) kl + ceekl] , (13a) 



D^^y{k) = Uf^^ ici2{k) + CQQ)kxky 



(13b) 



Dy^y (k) = 710^ [{ci2 (k) + 2c66) kl + ceekl] . (13c) 

The long-range Coulomb interaction renders the elastic 
coefficient Ci2 (but not the shear modulus cqq) nonlocal, 
so that Ci2 contains a diverging term 1/k. We shall 
write: 



Cl2 (k) = 



nk 



- Cl2, 



(14) 



where Ci2 is the weakly dispersive part of the elastic coef- 
ficient, and where the plasmonic (first) term on the rhs is 
due to the long-range nature of the Coulomb interaction. 

For the stripe state, Eq. ^ is no longer valid. In 
addition, all three elastic coefficients cii,ci2,C22 become 
nonlocal. We have in this case: 



(k) = riQ ^ [cii (/c) fc^ -f ceefc; 



^0 ^^y^ 



Dx.y (k 



rig ^ (ci2 [k] + cee) k^ky, 
Dy.y (k) = n-'^ [?22 {k) kl + ceekl] , 



(15a) 
(15b) 
(15c) 



where — — h Cy, with i,j = 1,2; and where 
we added to D^.x (k) a term Kky in order to take into 
account the bending rigidity of the stripes which, due to 
the small value of the shear modulus Cge in these systems, 
is quantitatively important over a sizeable region of the 
Brillouin zonei^. Using the fact that tiq = vjl-Kp'^ we 
finally obtain: 



no 



(16) 



We now want to discuss how one can evaluate the non- 
dispersive part Cij of the elastic coefficients. This will be 
the subject of the foUwing section. 



III. CALCULATION OF THE ELASTIC 
COEFFICIENTS IN THE HARTREE-FOCK 
APPROXIMATION 

In the Hartree-Fock approximation, a crystalline phase 
is described by the Fourier components {n (G)) of the 
average electronic density, where G is a reciprocal lat- 
tice vector. In the strong magnetic field limit where 
the Hilbert space is restricted to one Landau level, it 
is more convenient to work with the "guiding-center den- 
sity" (p (G)) which is related to (n (G)) by: 

(n(G)) =7V^Fjv(G)(p(G)), (17) 

where iV^ is the Landau-level degeneracy and 



Fn{G) 



74rO 



G2 



(18) 



is the form factor of an electron in Landau level N 
{L^jq {x) being a generalized Laguerre polynomial) . The 
magnetic field B = i?z is perpendicular to the 2DEG. 

The Hartree-Fock energy per electron in the partially 
filled Landau level is given bjs^ii^: 



E 



G 



(19) 

where the ^g,o term in this equation accounts for the 
neutralizing background of the ionized donors. The pa- 
rameter ly = Ne/Nip is the filling factor of the partially 
filled level, and we take all filled levels below TV to be 
inert. The Hartree and Fock interactions in Landau level 
N are defined by: 



i^(q) 
^(q) = 



1 -,^t^ 



dx ( 



(20a) 
(20b) 



X [L% (x^)]^ Jo {V2xq£y 



where Jq (x) is the Bessel function of the first kind. 

To compute the (p(G))'s, we first write this quantity 
in second quantization and in the Landau gauge A — 
{0,Bx,0) as: 



(P(G)> 



cjv x'^N,X-GyP 



(21) 

The average values (p(G)) are obtained by computing 
the single-particle Green's function (here and in what 
follows, Tt denotes the time ordering operator): 

G (X, X'., t) = - (TrCN,x (r) 4 X' (0) 



(22) 

whose Fourier transform we define as: 

^_ . G.(x+x')^^_^,_ G (X, X', r) , 



x,x 



(23) 
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so that: 



(p(G)) =G(G,r = 0- 



(24) 



We use an iterative scheme to solve numericahy^ the 
Hartree-Fock equation of motion for G (G,r) . For the 
undeformed lattice, we use the basis vectors: 



R-i = aoT] sin (93) x - 
R2 = aoy, 



aoricos ((p) y, 



(25a) 
(25b) 



where rj is the aspect ratio and ip is the angle between 
the two basis vectors. For the triangular lattice, rj — 1 
and ip = 7r/3. If we apply an elastic deformation u (r) 
to the lattice, the new lattice vectors are given by R' = 
nRi + mR2 + u (r) (where n, m are integers). We can 
write this expression as R' = nR[ 
the new basis vectors as: 



toRJ, if we define 



R'l 
R9 



aW sin [cp') 



aW cos{ip')y, 



(26a) 
(26b) 



The parameters Oq, 77' and ip' are functions of the original 
lattice and of the type of deformation considered. The 
reciprocal lattice vectors of the deformed lattice are eas- 
ily computed once these parameters are known. Then, 
the cohesive energy E{uq) of the deformed lattice can be 
calculated using the deformed reciprocal lattice vectors 
and Eq. ([9]). Under these circumstances, we find that 
the deformation energy per electron is given by: 



/ 



E{uo) E{uo = 0) 



(27) 



To find the elastic coefficients for the Wigner and stripe 
crystals, we need to consider the following deformations 
(note that the magnetic field and the number of electrons 
are kept fixedi^iH) : 

(i) A shear deformation with Ux (r) = uqi/ and Uy (r) = 
0: the strain tensors in this case are given by ex,x (r) = 
^v,v i'"^) — 0, and e^^y (r) = uo/2. The area of the system, 
5, is not changed by this deformation and the elastic 



energy is given by Fghear 
that the shear modulus cge is given by: 



j2 f 

1. ^ J shear 

C66 = lim no — — — , 
UQ^o duf. 



It then follows 



(28) 



where / = F/Ne is the deformation energy per electron. 
The parameters of the distorted lattice for this shear de- 
formation are given by: 

a'o = ao, (29a) 

ry' = T^y^l -I- Mo sin {(p) cos (tp) -I- Uq sin^ (ip), (29b) 

sin {(p) 



sin {(p') 



yjl + Uq sin (ip) cos {(p) + Uq sin^ (tp) 



(29c) 



(a) A one-dimensional dilatation along x, with 
Ux (r) = uqx and Uy (r) = 0: here, the strain tensors 



^x.x (r) = Uq, Cy^y (r) = 0, and ex,y (r) = 0, and the new 
area of the system is S" = (1 -|- Mq) 5 and F^x = ^SchUq. 
It then follows that the compression constant cn is given 
by: 



cii = lim no—-^ 



dx 



(30) 



while the parameters of the deformed lattice are given 
by: 



77' = T^y^l + (2mo + Mo)sin^ (V?), 

. , (1 -I- uo) sin ((p) 

sm yp ) — 



1 + (2uo + ul) sin^ [ip) 



(31a) 
(31b) 

(31c) 



The surface of the deformed lattice is S' = oIqT]' sin ((^') = 
S" (1 + Uq), so that the filling factor is now given by v' = 
i//(l + uo). 

(in) A one-dimensional dilatation along y with 
Ux (r) — and Uy (r) = u^y: now the strain tensors 
^x,x (r) — 0, ey,y (r) = mq and ex,y (r) = 0. The new area 
of the system is S" = (1 -I- mq) "S* and F^y = \Sc22'u^- The 
compression constant C22 is therefore given by: 



C22 = lim no- , 2 

Mo^o dug 



(32) 



On the other hand, the parameters of the deformed lat- 
tice are given by: 



ao = (l + uo)ao, (33a) 
-^J^ + {2uo + ul)cos^ (ip), (33b) 



V = 



(1 + wo) 



cos {(p') 



(1 + lio) cos {ip) 



+ {2uo + Uq) COS^ {ip) 



(33c) 



The surface of the deformed lattice is S' = a'Qf]' sin {p') = 
S [1 + Uq), so that the filling factor v' = v/ {1 + uq) . 

(iv) A two-dimensional dilatation with Ux (r) = uqx 
and Uy (r) = uoy: now, the strain tensors ex,x (r) = 
^y,y i^) ~ ""0 and Cx.y (r) =0. The new area of the system 
is S' = {1 + uof S and Fdxy = (cn -t- 2ci2 -I- C22) Uq. 
It follows that the combination cn + 2ci2 + C22 is given 
by: 



d fdxy 



Cn + 2ci2 + C22 = lim^ no — 



dui 



(34) 



For this case, there is no need to actually compute the 
energy of the deformed lattice since we can extract cn + 
2ci2 -l- C22 from the Hartree-Fock energy E/N^ given in 
Eq. (jl9p in the following manner. The area per electron, 
s, in the deformed lattice is s = (1 + uo)^ sq, so that {sq 
here is the area per electron of the undeformed lattice): 



Cn -I- 2ci2 + C22 = 4so 



^ fdxy 

ds^ 



(35) 
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The change in s causes a change in the fiUing factor, 
which is now given by: 



(36) 



Writing the HF energy as E/N^ = (fj) ^('^); have 
the relation: 



Cll + 2ci2 + C22 = 



e 



dv 



(37) 

Note that the long-wavelength Coulomb term '^""'^ 
must be added to Cii,ci2 and C22 that we compute in 
order to get Cii,ci2 and C22- 

Fig. 1 shows the expected quadratic behavior of the 
deformation energy as a function of uq, Eq. (j27|) . for 
a shear deformation in the small uq limit. In the one 
and two-dimensional compressions {i)-{iv), however, Eq. 
([27| leads to the addition of a non physical linear term 
in the dependence of the energies fdx 5 fdy 7 fdxy on ^0 
can be seen in Fig. 2. In the absence of deformation, the 
average electronic density is equal to that of the positive 
background. This neutrality removes the divergence of 
H (G) at G = in the Hartree-Fock energy of Eq. ([11]). 
When the electron lattice is dilated (but not the positive 
background), the electronic density no longer matches 
the density of the positive background and there is a 
restoring force that arises from this density imbalance. It 
is easy to show, assuming a density of the form of Eq. ([8]), 
that no linear term in uq arises when the interaction with 
the positive background is properly taken into account, 
and that the interaction with the background does not 
give rise to higher order terms in uq- Our Hartree-Fock 
procedure requires that the electronic and background 
densities be the same even for uq 0, which has the 
immediate consequence that we cannot directly compute 
the deformation energy using Eq. ((27l) . For all but the 
shear deformation, it is thus necessary for us to substract 
the linear term in Eq. (P7)l and to add by hand the 
long-wavelength Coulomb contribution of Eq. ^TU\i in 
order to get the correct elastic constants. Fig. 2 shows 
that a quadratic behavior for the deformation energy is 
recovered when the linear term is substracted. Note that 
the deformation energy in Fig. 2 does not contain the 
long- wavelength Coulomb contribution of Eq. (fTO)) . so 
that it can be either positive or negative. Our definitions 
in Eqs. ([501 [5^ of the elastic coefhcients are not 
affected by this procedure of removing the linear term in 
Mo since they involve the second derivative of the energy 
with respect to uq. 

It is instructive at this point to note that, for a triangu- 
lar Wigner crystal of classical electrons, the calculation of 
Bonsall and Maradudinii gives the following expression 
for the quantity A{u): 



2.5E-06 



2.OE-O61 




O.OE+00 



-0.01 



0.01 



FIG. 1: Deformation energy as a function of uq for a shear 
deformation: u (r) = mqi/x in a triangular Wigner crystal in 
Landau level A*" = with filling factor 1/ = 0.15. The square 
symbols are the HFA result while the solid line is a polynomial 
fit of order 2 (the linear term is negligible) . 



O.OE+00 



0.001 - 



0.000 

o 



-0.001 - 




8.0E-06 



-0.010 -0.005 0.000 0.005 O.OlU 

U„ 



FIG. 2: Deformation energy as a function of uo for a one- 
dimensional dilatation of a triangular Wigner crystal in Lan- 
dau level A'' = with filling factor u — 0.15. The square 
symbols are the HFA result (left axis) while the solid line is a 
polynomial fit. The dashed line (right axis) is the deformation 
energy with the linear term removed. 



A{iy) = -0.782133Vi^, 



(38) 
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and for the elastic coefficients: 



with the definitions: 



= -0.10892.3/^ (^) 
C66 = 0.015 56.3/2 (^-^). 



(39a) 
(39b) 



If we use Eqs. (|37p (with the relation ([4])) and take cge 
as given by Eq. (|39b[) . we find Ci2 — — Tcge, which is 
consistent with Eq. (|39ap . 



IV. DYNAMICAL MATRIX FROM THE GRPA 

We now turn our attention to the calculation of the DM 
of electron crystals in the GRPA method. In the strong 
magnetic field limit where the Hilbert space is restricted 
to one Landau level only, the Hamiltonian of the system 
is given by: 

iJ = ^ ^ li, (-k) D^.fi (k) up (k) , (40) 

k aS 

where a, /? = x^y and k is a vector restricted to the first 
Brillouin zone of the crystal. If we define the Matsubara 
displacement Green's function by: 

G„,^ (k, r) = - (r,M„ (k, r) up (-k, 0)) , (41) 

we find, using fi^ (...) — [i/, (. . .)] and the commuta- 
tion relation [u2;(k), Mj,(k')] — i^2(5k,-k', that this Green's 
function is related to the dynamical matrix by: 



Ga,/3 (k,ir2„) 



fi [172+^2^^ (k) 



(42) 



^^^^D,,,(k) -M^_i?,,^(k) 
Dy^x (k) Dx,x (k) 



where fin — 2T:n/T is a bosonic Matsubara frequency 
and 



(k) = -v/det [D{k)] 



(43) 



is the magnetophonon dispersion relation. 

We now define the following density Green's function: 

x[?^g, (k, r) = -N^ {Trpik + G, r) p (-k - G', 0)) , 

(44) 

where p = p — {p) . In the Generalised Random-Phase 
Approximation (GRPA), this Green's function is found 
by solving the set of equations^: 

[^^nSa,G' - A/g,G" (k)] Xg''!g' (k, «^n) = SG,G'(q), 



G' 



(45) 



Mc,G'(k)--2.(|;^)(p(G-G')) 

(k + G) X (k + G')' 



X sm 



X [H{G -G') - X {G -G') 
-i7(k+G') + X(k + G')], 



(46) 



and: 



Bg.g' (k) = 2z (p(G-G')) 

(q+G) X {ci+G')P 



(47) 



Diagonalizing the matrix M (k) and making the analytic 



continuation iQn 
the form: 



iS, we can write Xg'g' (^i ^) i^i 



E 



w. 



G,G 



(k) 



uj + i6 ~ LOi (k) 



(48) 



At small k, the pole oji (k) with the biggest weight 
Wq^'q (k) gives the GRPA magnetophon mode. We de- 
fine this pole as ujgrpa (k) and the corresponding weigth 
as (k). 

We now relate the displacement Green's function to 
the density Green's function using Eq. This last 

equation, coupled to Eq. gives the following relation 
between the density and displacement response functions 
(here F{k.) is the function defined in Eq. ^5]) . where for 
simplicity we now drop the Landau level index N): 



&^c^'(k,-) 



h{k + G)h{k + G') 
' F{k + G)F{k + G') 



(49) 



X ^(fca-l-G„)G„^^(k,w)(fc^ + Gy 



In deriving Eq. (|49|) . we have assumed that q • u (R) << 
1, so that a density fiuctuation can be linearly related 
to the displacement u (q) by Eq. ©. This is equiv- 
alent to assuming that the crystal can be described in 
the harmonic approximation so that only a knowledge of 
the dynamical matrix is necessary. To get Eq. (|49p. we 
have also assumed that h{r) — h (— r) so that h (q) is 
real. We can now use Eq. (|42p and the symmetry rela- 
tion Da,p (k) = Dp^a (k) to relate the density response 
function to the dynamical matrix: 



Ti{k) + i^T2 (k) 



[{L0 + t6f-u;l^{k) 

h{k + G)h(k + G') 
F(k + G)F(k + G')' 



(50) 
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where we defined: 

Ti (k) = -2 • [(k + G) X Z? (k) X (k + G')] • 2, (51a) 
r2(k) =2- [(k + G) X (k + G')]. (51b) 

For uj close to the magnetophonon resonance, we can 
write: 



Z(k) 



/i(k + G)/i(k + G') 



h uj + i5- u,np (k) F(k + G)i^(k + G') 

(52) 

where we defined the quantity: 



Z(k) 



ri(k) , .fir2(k) 



2cj,„p (k) 



2^2 



(53) 



Then, equating Eq. ((52)) with Eq. (l48l) for w close to 
<^GRPA (k), we obtain: 



Z(k) 



h{k + G)h{k + G') 



iS 



(k)F(k + G)F(k + G') 



(k) 



?'<5 - UJGRPA (k) ■ 



(54) 



Because a;,„p (k) must be equal to wgrpa (k) , we can 
finally write 



H ' F(k + G)F(k + G') 



(k), (55) 



or, taking the real and imaginary parts of this equation 
(we remind the reader that both functions /i(k) and F(k) 
are real): 



3fJ 



(GR.PA) 



G,G' 



w. 



{GB.PA) 



G.G' 



(k) 



(k) 



/i (k + G) /i (k + G') 



F(k + G)F(k + G') 
Ti (k)£4 ^ 

flWrnp (k) ' 

/i (k + G) /i (k + G') 



(k + G) (k + G') 
xEa (k) e. 



(56a) 



(56b) 



We can get rid of the unknown form factors /i (k + G) if 
we work with the ratio of the imaginary and real parts 
of the weights. We thus define: 



G.G 



(k) 





\jrAGRPA) 


(k)^ 
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\^AGRPA) 


(k)^ 



(57) 



(k + G) X i:>(k) X (k + G') 

~ nw„,p (k) (k + G) X (k + G') ■ 

A careful examination shows that, because uJmp (k) is 
given by the determinant of the dynamical matrix D (k), 
the quantity Ti (k) /tOrnp (k) is unchanged if all the com- 
ponents of the dynamical matrix are multiplied by some 
constant. Eq. (|57p is thus indeterminate. To avoid this 



problem, we replace uj^p (k) by ujgrpa (k) in Eq. (|57p 
Our final result is thus: 



G.G' 



-P (k + G) X i:>(k) X (k + G') 

" fil^GRPA (k) (k + G) X (k + G') ■ 

(58) 

Because D^.y (k) = Dy^^ (k), we need to choose three 
pairs of vectors (G, G') to get the components of the 
dynamical matrix. To be valid, the dynamical matrix 
obtained in this way must satisfy the equation: 



^GRPA (k) = -Vdet [D (k)]. 



(59) 



Eq. (f59|) provides a check on the validity of our calcula- 
tion. 



V. NUMERICAL RESULTS FOR THE WIGNER 
CRYSTAL 

In this Section, we illustrate the application of our 
method by computing the Lame coefficients for the trian- 
gular Wigner crystal in Landau levels iV = and N = 2. 
Fig. [3] shows the first two shells of reciprocal lattice vec- 
tors of the triangular lattice with Gi = (0, 0) . We take 
the vectors G and G' in Eq. ([58| on these first two shells. 
Not all combinations of vectors satisfy Eq. ([55]) . By ex- 
perimentation, we found that with a combination of the 
form [(G, 0) , (G', 0) , (G, G')] with G, G' ^ 0, this equa- 
tion is satisfied in the irreducible Brillouin zone shown in 
Fig. 3 to better than 0.05% for k£ < 0.3. We wiU thus 
stick to this type of combination for the rest of this paper. 




FIG. 3: First Brillouin zone of the triangular lattice with 
the irreducible Brillouin zone shown as the dark area. The 
arrows represent the reciprocal lattice vectors on the second 
shell while (1) corresponds to the vector G = 0. 

In the small- wavevector limit, the dynamical matrix 
of the Wigner crystal with a triangular lattice structure 
is given by Eq. (fTS]) . Using Eq. ([58]) . we can extract 
the elastic coefficients by fitting D^.fs (k) along the path 
ky = where: 
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0.04 




kj (k,=0) 



FIG. 4: (Color online) Component Dxx(k) of the dynamical 
matrix along computed for differents couples of vectors 
(G, G'). The numbers in the legend refer to the numerotation 
of the reciprocal lattice vectors shown in Fig. 3. 



D^,^ (k) = 
D.,v (k) = 0, 



e 



27r 



27r 

Dy^y (k) = —ceeklf, 



(ci2 + 2c66) (60a) 

V 

(60b) 
(60c) 



where cyi and cge are expressed in units of j kI^ . 

Fig. [4] illustrates one limitation of our method: the 
GRPA dynamical matrix is very much dependent on 
the choice of the couple (G,G'). Different choices give 
the same dynamical matrices D^^p (k) only in the small 
wavector limit kl < 0.1 as shown in Fig. 4 and, in this 
limit, the dynamical matrix element D^ r^ (k) is almost 
entirely dominated by the long-range Coulomb term (the 
first term in Eq. (|60ap ). It follows that different choices 
of (G, G') lead to quite different values of the elastic 
coefficients c\i even though Eq. ([59)1 is satisfied. The 
coefficient cee obtained from Z)j^^j^(k), however, is not 
affected by the long-range Coulomb interaction and ap- 
pears to be independent of the choice of (G,G'). Note 
that the dynamical matrix given by Eq. (j58p does not 
have the correct transformation symmetries of the trian- 
gular lattice. In cases where DaS (k) is needed in all the 
Brillouin zone, it becomes necessary to compute I?q./3 (k) 
in the irreducible Brillouin zone and obtain Da.ji (k) in 
the rest of the Brillouin zone by symmetry. 

To give an idea of the variability of the numerical re- 
sults with (G, G'), we show in Fig. [5] (for N = 2) and Fig. 
ini (for iV = 0) the coefficients Ci2 and cge extracted from 



the GRPA dynamical matric of the triangular Wigner 
crystal for different couples of vectors (G,G'). These 
coefficients are compared with those computed using the 
HFA described in Sec. III. We show the HFA results by 
a fuU line in Figs. 5,6. For both = and = 2, 
we find that the Hartree-Fock results for the coefficient 
C66 are extremely well reproduced by the GRPA method 
and, as we said above, do not depend on the choice of 
(G, G'). This is what we expect since the GRPA is the 
linear response of the crystal about the HFA ground state 
so that the coefficients Cy obtained from the two methods 
should be roughly equal, taking into account the various 
approximations made in deriving the GRPA dynamical 
matrix. The coefficient cee is easy to obtain, in view 
of Eq. (|60cp . since it is given by a one-parameter fit of 
the Dy y {kx: ky = 0) curve. The elastic coefficient ci2 
(which is related to the bulk modulus) is, on the other 
hand, much more difficult to obtain from Eq. (j60ap . In- 
deed, this elastic coefficient turns out to be very sensi- 
tive to how accurately the long-wavelength limit vk^i of 
Dx.x (k) in Eq. (|60ap is obtained by the GRPA numeri- 
cal calculation. (We here note that the GRPA dynamical 
matrix does contain the long-range Coulomb interaction 
discussed in Sec. II. The latter does not have to be added 
by hand as was the case for the elastic coefficients com- 
puted in the HFA.) As we see in Figs. 5(b) and 6(b), ci2 
is also very sensitive to the choice of the vectors (G, G') 
with one particular choice (2,3) reproducing the HFA re- 
sults almost exactly. The other two choices give very 
different values for ci2. In the absence of any criteria to 
choose (G, G') a priori, we would say that the GRPA dy- 
namical matrix cannot be used to make quantitative pre- 
dictions. The qualitative behaviour of the GRPA elastic 
coefficient Ci2 is consistent with that of ci2 computed in 
the HFA. 



If we exclude the domain v > 0.19 where our numer- 
ical results become noisy, we find that the average of 
the GRPA results for the three couples of (G,G'), as 
shown in Figs. 5(b) and 6(b), are in very good agree- 
ment with the HE calculation. In the absence of any 
criteria to choose the best couple (G, G'), this averaging 
procedure must be used to get qualitatively and quanti- 
tatively reliable results for the GRPA dynamical matrix. 
For 1/ > 0.19, the crystal softens and the quantum fiuc- 
tuations in u are important. There is a transition's into 
a bubble state with 2 electrons per unit cell at approx- 
imately ly = 0.22. We do not expect the assumptions 
underlying our method to be valid in this region. 
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FIG. 5: (Color online) Elastic coefficients (a) cge and (b) 
ci2 of the triangular Wigner crystal for Landau level A'' — 
2 computed using the different approximations listed in the 
legend. For the GRPA, the coefficients are computed using 3 
differents couples of reciprocal lattice vectors. The numbers in 
the legend correspond to the numeration of the vectors given 
in Fig. 3. The empty circles give an average of the 3 GRPA 
results. For cgg the different symbols are superimposed. 



FIG. 6: (Color online) Elastic coefficients (a) cge and (b) 
Ci2 of the triangular Wigner crystal in Landau level N — 
computed in the different approximations indicated in the 
legend. For the GRPA, the coefficients are computed using 3 
differents couples of reciprocal lattice vectors. The numbers in 
the legend correspond to the numeration of the vectors given 
in Fig. 3. The empty circles give the average of the 3 GRPA 
results. For cge, the different symbols are superimposed. 
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V 


cii (xlO-2) 


C12 (xlO-^) 


C22 (xlO-2) 


C66 (xlO^^^) 


0.42 


7.13 


-2.43 


-0.73 


4.35 


0.43 


5.91 


-2.47 


-1.15 


6.31 


0.44 


4.95 


-2.49 


-1.64 


7.08 


0.45 


4.21 


-2.50 


-2.15 


7.10 


0.46 


3.75 


-2.51 


-2.65 


6.21 



TABLE I: Elastic coefficients Wq ^c^j in units of / kH for the 
stripe crystal at various filling factors and in Landau level 
iV = 2. 



GRPA. We show in Figs. 8-10 the elements Dxx,Dj:y, 
and Dyy computed at filling factor v = 0.42 (in Lan- 
dau level = 2) along different directions in k-space 
together with the corresponding DM in the HFA element 
obtained from Eqs. ((6T|) with the coefficients of Table [D 
Similar results are obtained at other filling factors. No- 
tice that the bending coefficient K does not contribute to 
any of these curves. For each curve, Eq. is perfectly 
satisfied and the coefficient cee, which can be extracted 
from the GRPA function Dyy (kx, ky = 0), is in excellent 
agreement with the HFA results given in Table HI 



VI. NUMERICAL RESULTS FOR THE STRIPE 
CRYSTAL 

For the stripe crystal, the dynamical matrix is given 
by: 

+ Kk^i^), (61a) 
Dx,y{'k) = {-i^j—k.j;ky£'^ + — {ci2 + ce6)kxkyf, 

(61b) 

Dy, (k) = ^^klf + ^ {c,,klf + c,,klf) , 

(61c) 

and the elastic coefficients evaluated in the HFAl^ for 
Landau level N = 2 are listed in Table HI 




kj {k,=0) 



7 9 



4 2 13 5 
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FIG. 7: The first four shells of reciprocal lattice vectors of the 
anisotropic stripe cristal. 

The first 4 shells of reciprocal lattice vectors of the 
stripe crystal are represented in Fig. 7. From Eq. ([55)) . 
the vectors (G, G') must not be parallel otherwise the 
denominator in this equation vanishes. This forces us 
to use G in the second shell and G' in the fourth shell 
of reciprocal lattice vectors to evaluate the DM in the 



FIG. 8: (Color online) Component D^^{k) of the GRPA and 
HFA dynamical matrices of the stripe crystal computed along 
the direction ky = for partial filling factor u — 0.42 in 
Landau level N — 2, computed using 2 differents couples of 
reciprocal lattice vectors. The numbers in the legend refer to 
numbering of the reciprocal lattice vectors in Fig. 7. 

For the GRPA, Figs. 8-10 show results for the couples 
(G, G') that produce the maximum and minimum val- 
ues of the DM element. In all but the Dyy(kx = 0,ky) 
case, the HFA curve lies between these two results. For 
Dyy{kx = 0, fcy), oue of the GRPA curves almost coin- 
cides with the HFA result for kyt < 0.15. This is reas- 
suring for the validity of the GRPA method, but it also 
makes it impossible for us to find what part of the differ- 
ence between the GRPA and HFA is numerical and what 
part is physical (i.e. due to anharmonicity for example). 
We remark that, in the range k£ < 0.15, the GRPA re- 
sults are not numerically very different from the small 
k£ expansion of the dynamical matrix given in Eqs. (|6ip 
with the HFA coefficients. To the credit of our GRPA 
method, we add that the evolution of the different Dij's 
with filling factor is consistent with that of the corre- 
sponding elements calculated in the HFA as shown in 
Fig. 11. 
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FIG. 9: (Color online) Component DxyCk.) of the GRPA and 
HFA dynamical matrices of the stripe crystal computed along 
the direction ky = kx for partial filling factor u — 0.42 in 
Landau level N = 2, computed using 2 differents couples of 
reciprocal lattice vectors. The numbers in the legend refer to 
numbering of the reciprocal lattice vectors in Fig. 7. 
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FIG. 10: (Color online) Component DyyOi) of the GRPA 
and HFA dynamical matrices of the stripe crystal computed 
along the direction fc^, = for partial filling factor u = 0.42 
in Landau level N — 2, computed using 2 differents couples 
of reciprocal lattice vectors. The numbers in the legend refer 
to numbering of the reciprocal lattice vectors in Fig. 7. 



We thus see that the evaluation of the elastic coeffi- 
cients other than cge from the GRPA results seems haz- 
ardous for the stripe crystal. The curvature of the func- 
tions Dxx, Dxy, and Dyy in Figs. 8,9,10 is proportional to 
cii, ci2+ceG and C22 respectively. It is clear that the elas- 
tic coefficients extracted from these Dij are much bigger 
than those obtained from the HFA (the curvature of the 
HFA function is barely visible in the figures). These co- 
efficients also show very strong variation with the choice 
of (G, G'). An averaging of the GRPA results for differ- 
ent couples (G, G') would give a result closer to the HFA 
but the improvement would not be as dramatic as in the 
triangular lattice case. In fact, in the case of Dyy, we find 
that averaging over different choices of reciprocal lattice 
vectors does not bring any improvement to the numerical 
results. 
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FIG. 11; (Color online) Component Dxy{k) of the GRPA and 
HFA dynamical matrices of the stripe crystal computed along 



the direction kx = 
level N = 2. 



ky for different filling factors in Landau 



Finally, we remark that our GRPA results for 
Dxx {kx = 0, fcy) are dominated by a strong ky behavior 
indicating that the bending term K is absolutely essen- 
tial in the elastic description of the stripe crystal in Eqs. 



VII. CONCLUSION 

In conclusion, in this paper we have shown that is 
it possible to derive an effective dynamical matrix for 
various crystal states of the 2DEG in a strong magnetic 
field by computing the density response function in the 
GRPA. We have compared the dynamical matrix ob- 
tained in this way with the one obtained from standard 
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elasticity theory with elastic coefficients computed in the 
HFA. Our comparison was done for crystals with very 
different elastic properties, namely a triangular Wigner 
crystal and stripe crystal. Our motivation for deriving 
a dynamical matrix using the GRPA response consists 
in the fact that the latter has the advantage of giving 
the dynamical matrix directly without having to com- 
pute the elastic coefficients separately. Our comparison 
with the Hartree-Fock results showed, however, that the 
GRPA method must be used with care because of the 
variability of the results with the choice of the couples 
(G, G'). The shear modulus cge computed in the GRPA 
agrees very well with the one computed from the HFA, 
but the values of the other elastic coefficients Cij which 
are affected by the long range Coulomb interaction de- 
pend very much on the choice of the couples (G,G'). 
In some cases, as for a triangular Wigner crystal, an 
averaging procedure over different couples (G, G') im- 
proves the numerical accuracy of the method. In the 



long wavelength fc-^ <C 1 limit, however, the GRPA dy- 
namical matrix is a good approximation, both qualita- 
tively and quantitatively, and gives reasonable estimates 
for the elastic constants of the electronic solids that are 
in agreement with the static Hartree-Fock calculations. 
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